#include <iostream>
#include <pqxx/pqxx> 
#include <fstream>
#include <sstream>
#include <math.h>
#include <ctime>
#include <cstring>
#include <iomanip>

#define SIZE 46
#define YEAR 2011

using namespace std;
using namespace pqxx;

void apply_savitzky_golay4(double[], double[]);
void apply_lagrange_interp(double[], double[]);
double lagrange_interp(int, int, double, int, double, int, double, int, double);
bool get_4nearest_neighbor(int, double[], int[]);
void get_px_ndvi_lswi_evi_series( int, int, int, double, double[], double[], double[], int* );
void get_px_ndvi_series( int, int, int, double, double[], int* );
int count_0(double[]);
void write_series_to_file(string, int, int, double*);
void get_logs_for_sample_data();
bool is_not_water(int, int);

int main(int argc, char const *argv[])
{
    time_t t1 = time(0);
	// get thresholds
	/* file is named "thresholds.txt and formatted as:
		.(int)			threshold of sur_refl_b03 in xiao's method 
		.(int) (int)	min/max number of cloud to select data
	*/
	ifstream threshold;
	threshold.open("thresholds.txt");
	  	int th_xiao_sur_refl_b03;
	  	int min_cloud, max_cloud;
	  	threshold >> th_xiao_sur_refl_b03 >> min_cloud >> max_cloud;
	threshold.close();

	// get range of coordinate to process
	/* file is named "process_range.txt" and formated as:
		.(int) (int) 	top left
		.(int) (int)	bottom right
	*/		
	ifstream range2process;
	range2process.open("process_range.txt");
	  	int top, bottom, left, right;
	  	range2process >> top >> left >> bottom >> right;
	range2process.close();


	/*output file is name respectively and formated like:
		.(int) (int) (double[46]) 		line sample index_day_001 index_day_009...
	*/
	// for each pixel
	// count how many cloud cover, get 1, 2,..
	// get ndvi_seriess, lswi_series, evi_series
	// interpolate missing data
	// save to file ndvi_5cloud.txt

	int total_cloud;
	
	double ndvi_series[SIZE];
	double lswi_series[SIZE];
	double evi_series[SIZE];

	double interp_series[SIZE];
	double smooth_series[SIZE];

	int count = 0;
	int cloud_count[max_cloud-min_cloud]; for (int i = min_cloud; i < max_cloud; ++i) cloud_count[i-min_cloud]=0;

	for (int sample = top; sample <= bottom ; sample++ )
	{
		for (int line = left; line <= right; line++)
		{
			if (is_not_water(line, sample)){
				get_px_ndvi_lswi_evi_series(YEAR, line, sample, th_xiao_sur_refl_b03, ndvi_series, lswi_series, evi_series, &total_cloud);
				if (count_0(ndvi_series) >= min_cloud & count_0(ndvi_series) < max_cloud){
					write_series_to_file("ndvi.txt", line, sample, ndvi_series);
					write_series_to_file("lswi.txt", line, sample, lswi_series);
					write_series_to_file("evi.txt", line, sample, evi_series);
					
					apply_lagrange_interp(ndvi_series, interp_series); 
					write_series_to_file("ndvi_interpolated.txt", line, sample, interp_series);
					apply_savitzky_golay4( interp_series, smooth_series);
					write_series_to_file("ndvi_smoothed.txt", line, sample, smooth_series);

					apply_lagrange_interp(lswi_series, interp_series); 
					write_series_to_file("lswi_interpolated.txt", line, sample, interp_series);
					apply_savitzky_golay4(interp_series, smooth_series);
					write_series_to_file("lswi_smoothed.txt", line, sample, smooth_series);
					
					apply_lagrange_interp(evi_series, interp_series); 
					write_series_to_file("evi_interpolated.txt", line, sample, interp_series);
					apply_savitzky_golay4( interp_series, smooth_series);
					write_series_to_file("evi_smoothed.txt", line, sample, smooth_series);
					
					cout << "(" << line << "," << sample << ") total_cloud=" << total_cloud << ";"<< endl;
					cloud_count[total_cloud-min_cloud]++;
					count++;
				}
			}
		}
	}
	time_t t2 = time(0);
	cout << endl
		<< "Total " << count << " pixels in range " << min_cloud << "-" << max_cloud << ";\n" << endl
		<< "Total cloud distributed : \n";
	for (int i = min_cloud; i < max_cloud; ++i) cout << cloud_count[i-min_cloud] << "\t"; cout << endl;
    cout << endl
		<< "Start @ " << t1 << ";\n"
		<< "Finish @ " << t2 << ";\n"
		<< "Total processing time:\t" << t2-t1 << ";\n\n";

	get_logs_for_sample_data();

	return 0;
}



void get_px_ndvi_series(int year, int line, int sample, double th_xiao_sur_refl_b03, double ndvi_series[], int *total_cloud)
{
   *total_cloud = 0;
   string sql;
   stringstream sqlstream;
      
   try{
      connection C("dbname=mymodisdb user=demo1usr password=demo1pwd hostaddr=127.0.0.1 port=5432");
      if ( !C.is_open() ) { return; }

      sqlstream.str("");
      sqlstream << "\
         SELECT pixels.px_x, pixels.px_y, pixels.date_tkn, pixels.ndvi, pixels.lswi, pixels.evi, pixels.sur_refl_b03 \
         FROM public.pixels \
         WHERE px_x = "<< line <<" and px_y = "<< sample <<" \
         ORDER BY pixels.date_tkn \
         LIMIT 46 OFFSET "<< 46*(year-2009);
      sql = sqlstream.str();
      nontransaction N(C);
      result R( N.exec( sql ));

      int date = 0;
      double ndvi;
      for (result::const_iterator c = R.begin(); c != R.end(); ++c)
      {
         ndvi = 0.0;
         if ( c[6].as<double>() < th_xiao_sur_refl_b03 ){
            ndvi = c[3].as<double>();
         } else {
            (*total_cloud)++;
         }
         ndvi_series[date] = ndvi;
         date++;
      }
      C.disconnect ();
   }catch (const std::exception &e){
      cerr << e.what() << std::endl;
      return;
   }
}

void get_px_ndvi_lswi_evi_series(int year, int line, int sample, double th_xiao_sur_refl_b03, double ndvi_series[], double lswi_series[], double evi_series[], int *total_cloud)
{
   *total_cloud = 0;
   string sql;
   stringstream sqlstream;
      
   try{
      connection C("dbname=mymodisdb user=demo1usr password=demo1pwd hostaddr=127.0.0.1 port=5432");
      if ( !C.is_open() ) { return; }

      sqlstream.str("");
      sqlstream << "\
         SELECT pixels.px_x, pixels.px_y, pixels.date_tkn, pixels.ndvi, pixels.lswi, pixels.evi, pixels.sur_refl_b03 \
         FROM public.pixels \
         WHERE px_x = "<< line <<" and px_y = "<< sample <<" \
         ORDER BY pixels.date_tkn \
         LIMIT " << SIZE << " OFFSET " << SIZE*(YEAR-2009);
      sql = sqlstream.str();
      nontransaction N(C);
      result R( N.exec( sql ));

      int date = 0;
      double ndvi; double lswi; double evi;
      for (result::const_iterator c = R.begin(); c != R.end(); ++c)
      {
         ndvi = 0.0; 
         if ( c[6].as<double>() < th_xiao_sur_refl_b03 ){
            ndvi = c[3].as<double>();
            lswi = c[4].as<double>();
            evi = c[5].as<double>();
         } else {
         	ndvi = 0; lswi = 0; evi = 0;
            (*total_cloud)++;
         }
         ndvi_series[date] = ndvi;
         lswi_series[date] = lswi;
         evi_series[date] = evi;
         date++;
      }
      C.disconnect ();
   }catch (const std::exception &e){
      cerr << e.what() << std::endl;
      return;
   }
}

bool is_not_water(int line, int sample)
{
	stringstream sqlstream;
	string sql;	  
	try{
	  connection C("dbname=mymodisdb user=demo1usr password=demo1pwd hostaddr=127.0.0.1 port=5432");

	  sqlstream.str("");
	  sqlstream << "\
	     SELECT points.line, points.sample, points.land_cover \
	     FROM public.points \
	     WHERE line = "<< line <<" and sample = "<< sample;
	  sql = sqlstream.str();
	  nontransaction N(C);
	  result R( N.exec( sql ));

	  for (result::const_iterator c = R.begin(); c != R.end(); ++c)
	  {
	     if ( c[2].as<int>() != 0 )
	     	return true;
	  }
	  C.disconnect ();
	}catch (const std::exception &e){
	  cerr << e.what() << std::endl;
	}
	return false;
}

void apply_savitzky_golay4(double series[], double sg_series[])
{
	for (int i = 0; i < 4; ++i)
		sg_series[i] = series[i];
	for (int i = 4; i < SIZE-3; ++i)
		sg_series[i] = 1.0 / 231 * ( -21 *series[i-4]+ 14*series[i-3] + 39*series[i-2] + 54*series[i-1] + 59*series[i] + 54*series[i+1] + 39*series[i+2] + 14*series[i+3] - 21*series[i+4]);
	for (int i = SIZE-4; i < SIZE; ++i)
		sg_series[i] = series[i];
}

void apply_lagrange_interp( double series[], double lagrange_interp_series[])
{
	for (int i = 0; i < SIZE; ++i)
	{
		lagrange_interp_series[i] = series[i];
		if ( series[i] == 0)
		{
			int neighbor[4];
			if ( get_4nearest_neighbor(i, series, neighbor) )
			{
				lagrange_interp_series[i] = lagrange_interp(i, neighbor[0], series[neighbor[0]], neighbor[1], series[neighbor[1]], neighbor[2], series[neighbor[2]], neighbor[3], series[neighbor[3]] );
			}
		}
	}
}

double lagrange_interp(int x, int x1, double y1, int x2, double y2, int x3, double y3, int x4, double y4)
{
	double l1 = 1.0 * (x-x2)/(x1-x2) * (x-x3)/(x1-x3) * (x-x4)/(x1-x4); 
	double l2 = 1.0 * (x-x1)/(x2-x1) * (x-x3)/(x2-x3) * (x-x4)/(x2-x4);
	double l3 = 1.0 * (x-x1)/(x3-x1) * (x-x2)/(x3-x2) * (x-x4)/(x3-x4);
	double l4 = 1.0 * (x-x1)/(x4-x1) * (x-x2)/(x4-x2) * (x-x3)/(x4-x3);
	return y1*l1 + y2*l2 + y3*l3 + y4*l4;
}

bool get_4nearest_neighbor(int x, double series[], int neighbor[])
{
	int n = 0;
	for (int i = 1; i < 4; ++i)
	{
		if ( (x-i) >= 0 && series[x-i] != 0)
		{
			neighbor[n] = x-i;
			n++; if (n >= 4) break;
		}
		if ( (x+i) < SIZE && series[x+i] != 0)
		{
			neighbor[n] = x+i;
			n++; if (n >= 4) break;
		}
	}
	if (n == 4) return true;
	return false;
}

int count_0(double series[])
{
	int count = 0;
	for (int i = 0; i < SIZE; ++i)
		if ( !series[i] )
			count++;
	return count;
}

void write_series_to_file(string file, int line, int sample, double *series)
{
	char filename[100]; strcpy(filename, file.c_str());
	ofstream ofs; ofs.open(filename, ios::app);
	ofs << line << "\t" << sample << "\t";
	for (int i = 0; i < SIZE; ++i)
		ofs << setiosflags(ios::fixed) << setprecision(8) << series[i] << "\t";
	ofs << endl;
	ofs.close();
}

void get_logs_for_sample_data()
{
	int line, sample;
	double data[SIZE]; double tmp;
	ifstream ndvi_file, ndvi_interp_file, lswi_file, lswi_interp_file, evi_file, evi_interp_file;

	ndvi_file.open("ndvi.txt");	 ndvi_file >> line >> sample;  for (int i = 0; i < SIZE; ++i) { ndvi_file >> tmp;  data[i]=tmp;  }	 ndvi_file.close();
	cout << line << " " << sample << "\tndvi \t"; for (int i = 0; i < SIZE; ++i) cout << data[i] << "\t"; cout << endl;

	ndvi_interp_file.open("ndvi_smoothed.txt"); ndvi_interp_file >> line >> sample; for (int i = 0; i < SIZE; ++i) { ndvi_interp_file >> tmp; data[i] = tmp; }	ndvi_interp_file.close();
	cout << line << " " << sample << "\tNDVI \t"; for (int i = 0; i < SIZE; ++i) cout << data[i] << "\t"; cout << endl;

	lswi_file.open("lswi.txt");	lswi_file >> line >> sample; for (int i = 0; i < SIZE; ++i) {lswi_file >> tmp; data[i] = tmp;}	lswi_file.close();
	cout << line << " " << sample << "\tlswi \t"; for (int i = 0; i < SIZE; ++i) cout << data[i] << "\t"; cout << endl;

	lswi_interp_file.open("lswi_smoothed.txt"); lswi_interp_file >> line >> sample; for (int i = 0; i < SIZE; ++i) {lswi_interp_file >> tmp; data[i] = tmp;}	lswi_interp_file.close();
	cout << line << " " << sample << "\tLSWI \t"; for (int i = 0; i < SIZE; ++i) cout << data[i] << "\t"; cout << endl;

	evi_file.open("evi.txt");	evi_file >> line >> sample; for (int i = 0; i < SIZE; ++i) {evi_file >> tmp; data[i] = tmp;}	evi_file.close();
	cout << line << " " << sample << "\tevi \t"; for (int i = 0; i < SIZE; ++i) cout << data[i] << "\t"; cout << endl;

	evi_interp_file.open("evi_smoothed.txt"); evi_interp_file >> line >> sample; for (int i = 0; i < SIZE; ++i) {evi_interp_file >> tmp; data[i] = tmp;}	evi_interp_file.close();
	cout << line << " " << sample << "\tEVI \t"; for (int i = 0; i < SIZE; ++i) cout << data[i] << "\t"; cout << endl;
}